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A two-dimensional dense fluid of frictional grains is shown to exhibit time-chaotic, spatially het¬ 
erogeneous flow in a range of stress values, a , chosen in the unstable region of s-shaped flow curves. 
Stress controlled simulations reveal a phase diagram with reentrant stationary flow for small and 
large stress a. In between no steady flow state can be reached, instead the system either jams or 
displays time dependent heterogeneous strain rates j(r,t). The results of simulations are in agree¬ 
ment with the stability analysis of a simple hydrodynamic model, coupling stress and microstructure 
which we tentatively associate with the frictional contact network. 

PACS numbers: 83.80.Fg,83.60.Rs,66.20.Cy 


Discontinuous shear thickening is a ubiquitous phe¬ 
nomenon, observed in many dense suspensions L [21. 
Simulations for non-Brownian particle suspensions 3|-|5| 
as well as for frictional granular media |6|, LZ| have high¬ 
lighted the particular role played by frictional particle 
interactions. Discontinuous shear thickening implies a 
region of shear stress which (at finite Reynolds number, 
Ref. (H) is not accessible to a homogeneous system in 
stationary state. What happens, if the system is forced 
into this regime by prescribing the stress at the bound¬ 
ary in the unstable region? One possibility is vorticity 
banding @, corresponding to bands with different stress 
values at the same shear rate. However there is no clear 
evidence for persistent vorticity banding in experiment so 
far. Furthermore objections have been raised as to the 
possibility of vorticity banding as a stationary state: The 
pressure - in contrast to the shear stress - has to be the 
same across the interface of the bands; otherwise particle 
migration is expected to occur and thereby destabilise 
the interface. If stationary states are not accessible to 
the system, we expect to observe time-depend ent, inho¬ 
mogeneous states, either oscillatory or chaotic [10]. 

In this paper we show that spatio-temporal chaos oc¬ 
curs in a two-dimensional system of frictional granular 
particles, subject to an applied stress which is chosen in 
the unstable region of the flow curve. We present results 
from simulations and formulate a hydrodynamic model to 
derive a phase diagram and identify the regions of param¬ 
eter space, where time-chaotic, inhomogeneous solutions 
are to be found. 

We simulate a two-dimensional system of TV soft, fric¬ 
tional particles in a square box of linear dimension L 
as detailed in 0. The particles all have the same mass 
to = 1, but are polydisperse in size with diameters 0.7, 
0.8, 0.9 and 1 in equal amounts. Normal and tangential 
forces, / (n) and /W, are modeled with linear spring- 
dashpots of unit strength for both, elastic as well as 
viscous contributions. Thereby units of time, length 
and mass have been fixed Em. Flow curves for other 
visco-elastic parameters are presented in the appendix. 
Coulomb friction is implemented with friction parame¬ 


ter ij = 2, corresponding to the high friction limit. We 
expect the same qualitative findings as presented in this 
Rapid Communication for all values of /i > 0 and refer to 
a systematic study of the /r-dependence in ref. p]. In the 
stress-controlled simulations, a boundary layer of parti¬ 
cles is frozen and the boundary at the top is moved with 
a force crLe x , whereas the bottom plate remains at rest. 
In shear direction we use periodic boundary conditions. 

Constitutive equation - Previous work EE3 has 
revealed discontinuous shear thickening for a range of 
packing fractions close to the jamming transition. The 
flow curves for frictional granular particles are well rep¬ 
resented by the following constitutive equation 

j((j) = aa 1 / 2 — ber + ccj 2 , ( 1 ) 


where the term —ber is due to the frictional interactions 
and dependence on the packing fraction, </>, is imple¬ 
mented with a = a{(f >) UM- This gives rise to the 
phenomenology of the van-der-Waals theory for a first- 
order phase transition: jamming first occurs at the crit¬ 
ical point (f) c = 0.795 which also marks the onset of hys¬ 
teresis. A finite yield stress first appears at 4> a = 0.8003, 
while the generalized viscosity, 77 = cx/7 2 , diverges only 
at (fr] = 0.819. A similar, but not identical, sequence 
of characteristic packing fractions has been proposed in 


13}. Most remarkably is the existence (due to the in¬ 


dependent term) of a regime of shear stress which is un¬ 
stable < 0, corresponding to s-shaped flow curves. In 
strain-controlled conditions such an s-shape leads to dis¬ 
continuous shear thickening. In the van-der-Waals theory 
it corresponds to the coexistence region. 

Simulations in the unstable regime - These s- 
shaped flow curves are indeed observed in the simula¬ 
tions (see Fig. [1]) however only as transients and only 
in rather small simulation cells (TV < 24000). In larger 
systems the s-shape is slowly eroded and vanishes com¬ 
pletely above a certain system size. Instead a regime of 
continuous shear thickening develops. Closer inspection 
of the simulations in this regime reveals that no sim¬ 
ple steady-state is reached. Rather, the system displays 
time-chaotic and spatially inhomogeneous behavior. An 
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FIG. 1. Evolution of flow curves < 7 ( 7 ) with system size N. For 
small systems a pronounced s-shape is visible in transients be¬ 
fore the system jams. In large systems no steady flow occurs; 
average over the time-dependent flow results in continuous 
shear thickening (stress controlled except for N = 80.000). 



IVideo ll Four snapshots of the local stress, revealing large 
scale, time-dependent structures; <j> = 0.8035. 


example is shown in Video [Q which is a sequence of 
4 snapshots of the local-stress field (movies are given 
in the supplementary material [ijij]). The system does 
not settle into a time-independent steady state on the 
timescale of the simulations. Instead one observes time- 
dependent large scale structures, e.g. shear bands which 
seem to propagate in the principal stress direction, al¬ 
ternating with approximately homogeneous states and 
random large scale structures. In Fig. []] we show the 
corresponding time-dependent strain rate. One clearly 
observes irregular time-dependence with intermittent os¬ 
cillatory periods. 

Hydrodynamic model - In order to understand 
these time-dependent solutions and locate the regions of 


parameter space, where they can occur, we now formulate 
a hydrodynamic model, determine its stationary states 
and analyse their stability. Our starting point is the mo¬ 
mentum conservation equation in the form d t v x = d y a xy . 
For simplicity we only consider a one-dimensional model, 
allowing for a velocity v x in the flow direction, dependent 
on y only. In addition we introduce a variable w(y, t) for 
the internal state or the microstructure of the fluid. Such 
a variable has been introduced for many complex fluids, 
such as liquid crystals, entangled polymer solutions and 
colloidal systems Here we associate it with the fric¬ 
tional contact network. In the simplest model we only 
consider a scalar variable, representing e.g. the number 
of frictional contacts [§, [ij| , but are aware that a ten- 
sorial quantity, such as the fabric tensor, might be more 
appropriate. We assume that the microstructure variable 
relaxes to a stationary state dtW = (w — w*)/t which 
should vanish in the absence of stress w*(a —I 0 ) —> 0 . 
Furthermore dynamic rearrangements occur only due to 
driving, so that r _1 oc 7 . So far the model is the same 
as considered by Nakanishi et al. 0 for dilatant fluids. 
However, the coupling of stress and microstructure is dif¬ 
ferent for the frictional grains under consideration. The 
frictional contacts reduce the flow and hence the veloc¬ 
ity gradient. Starting from the constitutive equation for 
the strain rate of frictionless grains 70 = acr 1 / 2 + ccr 2 , we 
take the strain rate of our system of particles with fric¬ 
tion to be 7 = 70 — w. This completes the definition of 
the hydrodynamic model 

< 9 t 7 = dyer 
7 = 70 - w 

d t w =-^(w - w*). (2) 

Here we have introduced a proportionality constant, r = 
r /j, which can be fitted, when comparing simulations to 
the predictions of the model. For the stability analysis, it 
is irrelevant. Higher-order diffusive terms may be added 
to the stress and/or the w-equation. We have checked 
that the inclusion of such terms does not change the sta¬ 
bility analysis. Other hydrodynamic models of granu¬ 
lar fluids include velocity fluctuations [H], e.g. granular 
temperature which, however cannot explain effects due 
to friction. 

The model allows for two stationary states: The first 
one corresponds to stationary flow and is explicitly given 
by 


7 = 7o -w*; w = w*- cr = a 0 (3) 

Given that w* should vanish for vanishing shear, we take 
it as w* = ba , so that we recover the constitutive rela¬ 
tion for frictional grains, Eq. ©■ The second stationary 
solution accounts for the jammed state and reads 

7 = 0 ; w = 70 ; u = <tq 


( 4 ) 
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FIG. 2. Spatially averaged strain rate vs. time from nu¬ 
merically integrating the hydrodynamical model for differ¬ 
ent imposed stress. In the unstable region at intermedi¬ 
ate stress-values we observe oscillations and chaotic solutions 
(r= 10' 3 , 0 = 0.7975). 


For both stationary states, the stress a = a 0 is homo¬ 
geneous over the sample, in agreement with the Navier- 
Stokes equation which require a homogeneous stress in 
two dimensions and hence do not allow vorticity band¬ 
ing. 

Stability analysis - In order to study the stabil¬ 
ity of the stationary states, we consider small deviations 
5a, 5w ~ e nt e lky and linearise Eqj2] in 5a,5w. As ex¬ 
pected the stationary flow is unstable for < 0. Below 
0 C , this does not occur and we find two stable modes: 
a hydrodynamic one, fli oc — k 2 corresponding to the 
conservation of momentum, and a nonhydrodynamic one 
oc — ^ —> 0, corresponding to the relaxation of the 
microstructure, whose relaxation time becomes infinite 
as 0 —>■ 0 C and a —> a c . Above (j> c , the model predicts 
two stable stationary flow solutions, inertial flow at small 
stress and plastic flow at large stress. In between a gap 
of unstable stress values occurs, such that no stationary 
homogeneous flow is possible in this range of stresses. A 
typical flow curve in the range <f> c < (f> < (f> a is shown in 
the inset of Fig. Has line 1, indicating the unstable regime 
as red. The jammed state is only stable for 0 > 0 CT in 
the region where the constitutive relation yields a nega¬ 
tive j. A typical flow curve in the range (j> a < (f> < (f> v is 
line 2 in the inset of Fig. [3] with the stable jammed state 
marked in blue. For 0 > 0^, only stationary plastic flow 
and stationary jamming are predicted by linear stability 
analysis. The resulting phase diagram of the model is 
shown in the main panel of Fig. [3] 

Above the critical packing fraction, 0 C , a finite range 
of inaccessible stress values with ^ < 0 exists and gives 
rise to a corresponding range of unstable wavenumbers 
k 2 < k 2 = . which shrinks to 0 as 0 —> 0 C . Hence 


the model predicts an approximately harmonically oscil¬ 
lating state at the onset of instability, while well inside 
the unstable region more and more wavenumbers are un¬ 
stable so that one expects a broad range of frequencies to 
be present in the spectrum. These expectations are born 



FIG. 3. Inset: Typical stationary states (flow curves) in the 
range 0 CT > 0 > 0 C (1) and 0 > 0 CT (2). Unstable regions 
are highlighted in red (time dependent flow) and blue (jam¬ 
ming). Main: Phase diagram following from the linear sta¬ 
bility analysis; the two generic flow curves, displayed in the 
inset, correspond to the paths, denoted by 1 and 2. 


out by numerical integration of the partial differential 
equations in order to obtain the full nonlinear dynamical 
evolution. Close to 0 C , the oscillations are approximately 
harmonic, while we find oscillating and seemingly chaotic 
solutions at larger packing fractions (see Fig. [2]). 

Comparison with simulations - To check the 
predictions of the above analysis, we performed stress 
controlled simulations (for technical details see Ref. 0 ) 
along the paths 1 and 2 in the phase diagram. The time- 
dependent strain rate along path 1 is shown in Fig. 0}), 
for three different stress values. The lowest one shows 
stationary flow in the Bagnold regime, the chaotic time- 
dependence with oscillatory components is represented 
by the two red curves for intermediate stress and larger 
stress gives rise to stationary plastic flow (not shown in 
Fig. Hi). Similarly in Fig. H J we show the strain rate as 
a function of time for path 2. In addition to the steady- 
state flow, and the oscillatory flow, there are stationary 
jammed states, where the initial flow ceases after a cer¬ 
tain amount of time. 

To get a quantitative measure of the irregularity in 
the time dependence of the states, we have computed 
the power spectrum C(w) as the Fourier transform of 
the strain rate auto-correlation function (Fig. [5]). In the 
stable regimes, the power spectrum is C ~ u; -2 , sug¬ 
gesting simple exponential correlations and linear noise. 
In the unstable regime this background spectrum is su¬ 
perposed by additional and irregular complex structures. 
This is a strong indication of truly nonlinear chaotic dy¬ 
namics. Noteworthy is also the strong peak at higher 
stresses. This corresponds to the fast oscillations visible 
in Fig. Hi. 

Finite-size effects - Snapshots from video |T] indicate 
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FIG. 4. Strain rate 7 (t) as a function of time for several values 
of stress, corresponding to path 1 (top, <j> = 0.7975) and path 
2 (bottom, (f> = 0.8035) in Fig. 0 the sudden drops of the 
strain rate in the time dependent flow curve indicate that the 
system nearly gets jammed. 


the build-up of large-scale coherent structures. A large 
correlation length has also been observed in Ref. [3| in the 
context of continuous shear thickening in non-Brownian 
particle suspensions. These correlations are also in line 
with the strong finite-size effects observed in the flow 
curves of Fig. [TJ Indeed, lowering the system size the 
oscillatory state acquires a finite lifetime and the system 
jams. This is because in small systems there is a finite 
probability that large strain rate fluctuations towards 
7 —> 0 lead into a stable jammed state (see Fig. SJo for an 
example of such an excursion). Similarly, the appearance 
of s-shaped flow curves is due to system-size dependent 
strong strain rate-fluctuations towards the jammed state 
with 7 = 0 . 


Conclusion - We discuss stress-controlled driving 
of a granular system that undergoes discontinuous shear 
thickening. In particular a regime is identified where the 
system does not settle into a time-independent steady- 
state. Instead, it displays spatio-temporal oscillations 
and chaotic behavior as a novel possibility to adopt to 
stress in the unstable parts of the flow curve. Recent 
experiments on corn-starch reveal very similar unsteady 
flow where theory predicts discontinuous shear thicken¬ 
ing [l2|. 


Simulations reveal a phase diagram that has a char¬ 
acteristic re-entrant form with steady-state flow at small 
and large stresses. At intermediate values of stress either 
time-dependent states are observed or the system settles 
into a non-flowing jammed state, depending on stress and 
packing fraction. For </> c < (f> < cf> a only time dependent 
solutions are observed, whereas for (f> a < <j> a sequence of 



FIG. 5. Power spectrum C(uj) of strain rate fluctuations at 
constant stress from low (bottom) to high (top) values of a 
(shifted vertically for clarity of presentation); in the stable 
region C ~ w -2 , whereas in the unstable region additional 
complex structures are superimposed on the u ;~ 2 decay. 


inertial flow, chaotic flow, jammed state and plastic flow 
are seen for increasing values of tx, until at <f> = (j)^ only a 
transition from the jammed state to plastic flow remains. 

We also present a hydrodynamical model, coupling 
stress to a microstructural observable. Within linear 
stability analysis we recover the detailed features of the 
phase diagram as obtained from simulations. In the un¬ 
stable region, the model predicts either oscillating or 
time-chaotic flow. 

In future work we plan to quantify the spatio-temporal 
correlations that are visible in the snapshots and com¬ 
pare them with length-scales determined in [cjj from ve¬ 
locity correlations. We furthermore aim to better under¬ 
stand the microstructural observable, check whether it 
can be associated with the contacts which are blocked by 
Coulomb friction and explore the possibility of a tensorial 
observable, such as the fabric tensor. 
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Appendix 

In this appendix we deal with flow curves of frictional 
grains at fixed packing fraction when the particles’ elastic 
and viscous damping constants change. 

In the main article we use unit strength for elastic and 
viscous contribution in the linear-spring dashpots (fc and 
77 , respectively). Units of time and energy dissipation 
when particles interact are set by this choice. In partic¬ 
ular, the coefficient of (normal) restitution, e^- n \ and the 
binary collision time, t, are set. Table [T] summarizes the 
normal coefficient of restitution and the binary collision 
time for the parameters that we discuss here. We tune 
both contributions, k = k = k^ and 77 = r 
independently. Both, and t, just serve as a guidance 



FIG. 6. Scaled flow curves for (j> = 0.80, viscous damping 
parameter 77 = 1 and varying elastic constant k. 


show a choice of packing fractions with phenomenology 
similar to Fig. [ 6 ] Lower <j> leads to flow curves similar to 


TABLE I. Normal restitution coefficient t ^ and binary col¬ 
lision time t for different k and r/. 


k 

V 

e (") 

t 

1 

1 

0.305 

2.375 

1 

1/2 

0.569 

2.26 

1 

1/10 

0.895 

2.22 

1/2 

1 

0.163 

3.63 

2 

1 

0.44 

1.4 

10 

1 

0.7 

0.7 


on properties of pairwise collisions and do not respect the 
large packing fraction. Also the important frictional con¬ 
tribution which implies tangential restitution is not 
characterized by these quantities. The normal restitu¬ 
tion can be computed easily while the tangential part is 
of rich and complicated behaviour due to its dependence 
on the impact velocities [ 20 ]. 

The following data shows flow curves of a system with 
N = 8000 particles. Fig. [ 6 ] shows scaled flow curves for 
fixed (j) = 0.80 and 77 = 1. When k decreases (e^") de¬ 
creases), the discontinuity shifts towards small strain rate 
until it is missed out by our simulation. Note that the 
numerical effort for low strain rate is much larger as for 
large strain rate. An increasing k (increasing e*™') shifts 
the discontinuity towards larger strain rate which makes 
the discontinuity smaller until it vanishes. At k = 10 
the transition from inertial to plastic flow is smooth and 
without shear thickening. The phenomenology of varying 
k seems in that range similar to a change of the packing 
fraction. The latter was studied in previous work Q. 
A systematic study of variations of all three parameters, 
(j), k and 77 , is out the scope of this article. In Fig. Q] we 



FIG. 7. Flow curves for varying packing fraction (j> across the 
transition with viscous damping parameter 77 = 1 and elastic 
constant k — 1. 

those with large and large 4> show the same behaviour 
as small e^ n \ The variation of the viscous damping pa¬ 
rameter 77 is shown in Fig. [5] The phenomenology is the 
same as above: decreasing e ^ (larger 77 ) leads to a shift 
towards small strain rates and an increasing coefficient 
of restitution (smaller 77 ) shifts the discontinuity towards 
larger strain rate until it vanishes. Then the transition 
between inertial and plastic flow is smooth and without 
shear thickening again. 
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FIG. 8. Flow curves for cj> — 0.80, elastic constant k = 1 and 
varying viscous damping parameter rj. 
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